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I . INTRODUCTION 



The problem of identifying the dynamics of a physical 
system from its input-output behavior is a fundamental issue 
in modeling and control. The most general approach is based 
on the assumption of a mathematical model, often linear and 
time invariant, and its parameters are determined by 
optimization of an error criterion. 

A wide class of identification algorithms exists in the 
literature. Several books have been written on this matter, 
see for example [Ref. 1, 2]. Of particular interest are the 
recursive algorithms, where the estimated parameters are 
updated on-line based on new input-output data collected at 
each sampling instant. In most of the approaches on recursive 
identification the parameter estimates are updated for every 
new set of data points; as a consequence, the estimate at any 
time is updated up to the current data. 

A different approach is based on Block Processing (BP) , 
where the parameters are recursively updated on the basis of 
the data within nonoverlapping intervals of time. Although 
this introduces a delay between the time we measure the input- 
output data and the time we use it in the estimation, it has 
the advantage that time averaging within the time interval 
reduces the effect of disturbances. 
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Block Processing techniques have been investigated by 
several authors. Of particular interest for this thesis is 
the approach developed by Mansour [Ref. 3] in the adaptive 
filtering contest, where the parameters are updated on the 
basis of the DFT (Discrete Fourier Transform) of the data 
within each time interval. 

In this thesis an adaptive identification algorithm 
operating on the DFT (Discrete Fourier Transform) of blocks of 
input and output data is introduced. In particular, the 
algorithm identifies the frequency response of the system's 
linear model by applying several recursive estimation 
techniques and adapting them to the frequency domain approach. 
We will investigate estimation methods based on the Projection 
Algorithm and Recursive Least Squares [Ref. 4] and the results 
will be compared. Even though Recursive Least Squares is a 
complex and time consuming method, it is preferable when 
accuracy and speed of convergence are needed. 

This thesis is organized as follows. Chapter II 
introduces the concept of Block Processing and develops two 
algorithms for recursive identification in the frequency 
domain. Proof of convergence is also part of this chapter. 
Simulation studies are introduced in Chapter III showing the 
effectiveness of the algorithms and the conclusion is in 
Chapter IV. 
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II. IDENTIFICATION OF LINEAR SYSTEMS 



A. STATEMENT OF THE PROBLEM 

Consider an LTI (Linear Time Invariant) system with input, 
output signals given by u(t) and y(t) respectively. The 
problem we address is the estimation of the parameters of a 
linear model in order to predict future values of the output 
y given past measurements of the input and output signals. 
The general block diagram is given in figure 2.1. The block 
z' 1 denotes time delay by one clock pulse. 




Figure 2.1 Block diagram of the system 



This figure shows the general scheme where y(t) is the 
prediction of y(t) based on the past input, output, and the 
estimated model. The purpose is to determine a model for the 
plant so that the error e(t) is miniminized. In particular, 
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the model is restricted to be an LTI system in either one of 
the two forms: 

ARMA (Auto Regressive Moving Average) 

y ( t) =a 1 y( t-l) + . . . +a n y( t-n) +i? 1 u( t-l) +. . .b n u( t-n) (2.1) 

where a { , are constants; or 

CONVOLUTION form, where the output input signals are related 
by the convolution sum. 



y(t) =h ( t) *u( t) 



( 2 . 2 ) 



y(t) =]P h( t) u(t-x) 

T = -oo 



( 2 . 3 ) 



The problem of adaptive frequency response identification 
mainly has been concentrated on the approximation of infinite 
impulse response systems by finite impulse response models. 
In this thesis, we will consider the identification of stable 
plants only, for which the impulse response h(t) decays to 
zero as time increases. For this reason, the infinite impulse 
response h(t) can be approximated by a finite sequence with 
length N. In this way, the moving average (MA) model of the 
following form is obtained. 
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y(t) =h( 0) u(t) (1) u(t-l) +-+h (N-l) u(t- (N-l) ) (2.4) 



where h(t), for 0 < t < N-l, is the truncated impulse 
response, We choose the block length N to be a power of 2 in 
order to apply FFT (Fast Fourier Transform) techniques. 

The approach we consider is based on a frequency response 
framework. In particular, the estimated model is computed 
from the Fourier Transform of blocks of data collected from 
measurements of the input and output. 

There are several reasons which motivate this approach. 
In many cases of interest, a model which is accurate on a 
desired frequency range only is identified for the plant. For 
example, if we want to model the low frequency spectrum of the 
system with the frequency response approach, we select weights 
in order to select frequencies of interest. 

In the next section, a recursive algorithm for the 
identif icaton of the frequency response is introduced. 

B. RECURSIVE FREQUENCY DOMAIN APPROACH 

Block Processing (BP) techniques for on line 
identification of linear models have been investigated by 
several authors [Ref. 3]. With this approach the estimated 
model parameters are sequentially updated on the basis of 
blocks of data, rather than at each data point. In this 
section, we investigate a BP adaptive algorithm which operates 
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on the Fast Fourier Transform (FFT) of blocks of input and 
output data. 




There are some advantages in dealing with data by blocks 
rather than individually. Block processing averages the data 
within the block, resulting in better noise rejecton. In 
order to illustrate the BP approach in the ARMA model case, 
consider the model in (2.1) and the data y k (t) , u k (t) within 
the k-th time block, i.e. 

y k { t) =y{kN+t) 
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u k ( t) =u ( kN+t ) 



(2.5) 



for t = 0, 1, ... , N-l. If the block size N is larger than 
the system order n, we can write (2.1) as 

y k ( t) =a 1 y Jt ( t-l) +a 2 y*( t-2) +b 1 u k ( t-l) + b 2 u k (, t-2) (2.6) 

for a second order case with 2 < t < N-l. It is easy to see 
that (2.6) can be written in convolution form as 

y k ( t) =a ( t) ®y k ( t) +b{ t) ®u k ( t) nztzN-1 (2.7) 

with a and b sequences in R n 

a=[0 a x a 2 - a n 0 - 0] 

£>=[0 b x b 2 - b n 0 - 0] (2.8) 

and ® denoting circular convolution. 

In the approach, the well known relationship between the 
DFT and the circular convolution in the frequency domain is 
exploited. By defining { Y k (l) , 1=0,1, ...,N-1 } and { U k (l) 

, 1=0,1, ...,N-1 } as the DFT of the output and input data 
y k (t) ,u k (t) in the k-th block, we can define ? k (l) as 

Y k (l)=A(l) Y k (l) +B(l)U k (l) 2=0,1, — , N-l . (2.9) 
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Clearly from (2.7) and (2.9), we can write the important 
relationship 

y k (C)=y k (C) n<.tzN-l (2.10) 



with y k (t) = IDFT [ Y k (l)], and IDFT denoting inverse DFT . 

A linear estimation algorithm for A(l) and B(l) in (2.9) 
can be determined by defining an error signal 



e k (t) 



y k (t) -y k ( t) 



0<. tzn-l 



( 2 . 11 ) 



where y k (t) = IDFT[ A k (l) Y k (l) + B k (l)U k (l) , 1=0, 1, . . . ,N-1 ] and 
A k , B k are estimates available at the k-th block. From the 
above definitions, we can write the recursive estimates for A, 
B by the algorithm. Let 

©,(J)=[4<j),e t (i)] r 



X x U) =[Y k {l) ,U k (l)] T (2.12) 

for 1=0, 1, ..., N-l , k = 0, 1, 2, ...; also let 

E k (l) =DFT[e k ( t)] (2.13) 

with e k as in (2.11). Then the recursion is such that 
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(2.14) 




For this algorithm we can show the following: 



(a) || 6 k+1 (l) - 0(1) II 2 < II 0 k (l) - 6(1) II 2 



(2 . 15) 



for all 1=0, 1, 2, N-l, and for all k; 

(b) If the input and output data u k (t), y k (t) are bounded for 
all k and t, then 



In this result we see that the parameter error decreases 
with time (eqn. (2.14)) and the error between the model and 
the plant e k (t) tends to zero, i.e. the model output tends to 
follow the plant output. 

Proof: Let us define the parameter error 



Then from (2.14) and (2.17) we can write the recursion 



e k (t)=0 for all 0<.tzN-l 



(2 . 16) 



& k (D =® k u) -e k (i) 



(2.17) 



& k . 1 (l)=Q k (l)- V i k (X t k (l)E k (l)) 



(2 . 18) 



where 






1 



1 + j (\U k ( 1 ) | 2 +\Y k ( 1 ) | 2 ) 



(2.19) 
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By taking the magnitude square of both sides of (2.18), we 
obtain 



= [&' T k (l) -\l k E ,T k (l)xZ(l)] [& k (l) -\i k X k (l)E k U)] 

=p k (i)\\ 2 -\i k E* k T (i)x£(i)& k (i)-\i k e k r (i)x' k (i)E k (i) 

+ \i 2 k E k T xZ(l)X'k(DE k (l) ( 2 . 20 ) 



By taking the sum of all frequency components in both sides of 
(2.20), we obtain the expression 

1 =0 1=0 1=0 



& , k T (DX k (l)E k (l) 

1=0 



N-l 



1=0 



\i k E k T (l)X k (l)X k (l)E k U) 



( 2 . 21 ) 



Applying Parseval ' s theorem, we can relate data in the 
frequency domain to the corresponding time domain as 



10 



( 2 . 22 ) 



N - 1 

E 

2=0 



N - 1 



[5; T (j>] t^ t r (i)e Jt (j>] «£ § k (ot k (t) 

t=o 



N - 1 



E 



tf-1 



e Jt (t)e^(O=5]|e Jt (t)| 2 ^0 

t-2 



where we define 



e ic ( t) =JDFr(zJ ( J ) J > ] 



Equation (2.23) comes from the fact that 



e^( t) =0 



Os t<.n-l 



e k (t) =y k (t) -? k (t) =e k (t) nztzN-1 



and 



c k (t)=e k (t) 



Still using Parseval's theorem, we can write 



EP*<-»r-pj 



j=o 



2 



(2.23) 



(2.24) 



(2.25) 



(2.26) 



(2.27) 



(2.28) 
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After some manipulations and the fact that | X k ( 1 ) | 2 < 1 for 
all k and 1 = 0, 1, N-l, we can bound the parameter error 
at the end of the (k+l)th block as 

IMI 2 * IM 2 -2NN + NN 2 (2 - 29) 



and therefore 




; -Nl g *ll 2 



(2.30) 



Finally, since || © k II > 0 for all k and is a nonincreasing 
sequence as in (2.30) , the increment ^ k || e k || 2 must tend to zero 
as k tends to infinity. Therefore 



lim 



lie J | 2 



k ~°° l+MAX§Y k (l)\ 2 +\U k (l)\ 2 ] 



= 0 



(2.31) 



which proves the result. 



The significance of this result is that we can estimate 
the parameters of the given system (in terms of A k (l), B k (l)) 
based on the frequency content of each block of data. 
Furthermore, the estimation is recursive, and the convergence 
of the prediction error e(t) to zero is guaranteed at least in 
the ideal case. It will be shown in the simulations that even 
in the presence of measurement noise, the convergence of the 
error to small values is still satisfactory. 
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In the implementation of the estimation algorithm 
introduced above, we need to compute the DFT of the "windowed" 
error term e k (t) in (2.25) and (2.26). We can see that in the 
implementation two operations of IDFT are required to compute 
y k (t) . Using DFT again, E k (l) is obtained, which is used in 
the recursion (2.14). 

An algorithm, which does not require this sequence of IDFT 
and DFT in its implemention, is introduced next. The new 
algorithm also offers the advantage of being able to use 
different recursive estimation techniques (such as the 

Recursive Least Squares) which exhibits faster convergence. 
However, this is obtained at the expense of added complexity. 

In order to introduce the technique, let us recall the 
signal 

y k (t) =a{t)®y k { t) +b{t)®u k {t) (2.32) 

with a, b being the plant parameters, and the equality 

y k (t) =y k (t) nztzN-i (2.33) 

Going to vector notation, let us write (2.33) in the following 
form 



13 



(2.34) 



0 






' y*(o) 


y k (n) 




[0] [0]' 


9 k in) 






.[0] [I], 


: 


y k (N-2) 






y k (N- 2) 


y k (N- 1) 






? k (N-l) 



=h9 k 



where [0] and [I] denote blocks of zeros and the identity 
matrix respectively. Definition of the Fourier matrix F as 



F =e 

r m, n c 



_ • / 2 nmn , 



(2.35) 



allows the computation of the DFT of a sequence as a matrix 
operation. In this way, we obtain from (2.34) 

Fy k =Fhy k (2.36) 

and therefore, after simple manipulations, the following 
equation, 

Fy k ={FhF~ x ) (Fy k ) (2.37) 

It is easy to see that the vectors Fy k and F y k are the arrays 
of DFT coefficients of the respective sequences y k and y k , 
which yields 

Y k =HY k (2.38) 



Recall from (2.32) that 
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Y k (l) =A k ( 1) Y k (l) +B k (1) U k (l) 



(2.39) 



for 1=0, 1, 2, N-l. Let us define YU as the following 

expression 



YU= 



Y k { 0) 0 

0 Yj. ( 1 ) 



0 



0 

0 

0 Y.(N- 1) 



5 u k { 0 ) 0 

i 0 U k { 1) 

: : 0 



0 

0 

0 UAN- 1) 



(2.40) 



and combine (2.38) and (2.39) to yield 



**= (■ HYU) 



■ A k ( 0 ) 
A k ( N- 1 ) 
B k (0) 
BAN- 1) 



(2.41) 



The recursive estimation algorithm is derived from (2.41) . In 
particular, we know that (2.41) can be broken into a sequence 
of linear equations of the form 

Y k U) =<J>J(i)0 1=0,1, 2, . . . ,N-1 (2.42) 



where 
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0=[A*(O) - A k (N- 1) £*(0) - B k (N-l) ] T (2.43) 

is the vector of unknown parameters, and 

<&Ul) = [H{1,0) - H(1,N-1) ] (YU) 

= [H(1,0) Y k (0) - H{1,N-1) Y k {N-l) 

,H(l,0)U k (0) - H(l,N-l)U k (N-l)] 

k k (2.44) 



In (2.44) the coefficient H(i,j) are constant (in the sense 
that they do not depend on the block k) and are determined by 
the window matrix h only. 

From (2.42) we can use any recursive algorithm to estimate 
0. In particular, a Recursive Least Squares [Ref. 2] 
algorithm applied to complex data is going to yield the 
following recursion 






Pt* k (i) (y je (i)-oj(i)0j) 

i+$ T k (i)p£$ k (i) 



(2.45) 
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(2.46) 



k k 1 +Ql(l)pfa k (l) 

for 1 = 0, 1, 2, . .., N-l, with initial conditions at the 

beginning of each time block given by 

&° k =6 N k _i ; Pk=° 2 I (2.47) 

with I the identity matrix and a 2 an arbitrary constant 
parameter. In general, the constant a is chosen to be a large 
value, in comparison with the order of magnitude of the 
parameters 0 . 

Although far more complicated to implement, the algorithm 
(2.45), (2.46), and (2.47) has better convergence properties 

than the one shown previously. 

In the next section the two algorithms are applied to the 
identification of a linear time invariant system. 
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SIMULATION STUDIES 



A. INTRODUCTION 



The 


estimation 


techniques 


discussed 


in 


the previous 


chapter 


have been 


simulated, 


and used 


to 


estimate the 



frequency response of a system. In addition, two computer 
programs written in MATLAB have been developed. The first one 
is the main program given in Appendix A which simulates the 
entire system in an iterative fashion. The main program 
estimates the value of the frequency response of a given 
system at a particular frequency and compares it with the 
corresponding value of the original system. 

The second program is a function type of subroutine which 
implements the Recursive Least Squares algorithm [Ref. 2]. 

The proposed identification method has been tested using 
several different inputs in the program in order to simulate 
and analyze the effect of various excitaton conditions. In 
the next section, a first order recursive difference equation 
is used as an example. 

B. SIMPLE DERIVATION 

The system being used as an example is a first order 
system described by the difference equation. 
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y(t) =0 . 3y ( t-1) +5 . Ou ( t-1) 



(3.1) 



In like manner (3.1) can also be written in terms of the 
impulse response, 

oo 

y(t) =h( 0) u ( 0) +h( 1) u(t-l) + . . . = h(m) u( t-m) (3*2) 



Since the system is stable, its impulse response decays to 
zero, and (3.2) can be approximated by a finite sum 

y(t)=h(0)u(t) +h (1) u ( t-1) +. . . +h (N-l) u ( t-N+1) (3.3) 



In the followig simulations N is chosen as N = 32 data points. 
The transfer function H(z) of the original system 



* (z)= Zi4 = _L_ 

U(z) z- 0.3 



(3.4) 



yields its frequency response 

H(e^)=— i (3.5) 

e j6 -0 . 3 

Due to the use of FFT methods, we estimate the frequency 
response at discrete frequency points 0 = 27rl/N, 1 = 0, 1, 2, 
..., N-l. For this purpose, let us define 
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■ 4-00 



(3.6) 



H(l) = h (n) e 




and its finite approximation 



H(l) =J^h(n)e 



N - 1 




(3.7) 



The next step is to calculate the frequeccy response of the 
original system. 

C. STRUCTURE OF THE PROGRAM 

The algorithms introduced in the last chapter have been 
tested in the identification of a discrete time system. In 
particular we look at the problem of identifying the frequency 
response of the system by fitting a finite impulse response 
system to the input-output data. 

The software developed consists of two programs: a main 
program which produces the input output data used in 
identification, and a subroutine which implements the 
Recursive Least Squares identification shown in Chapter II. 

The effectiveness of the estimation is assessed on the 
basis of two criteria: the error between actual output of the 
system and predicted output, and the error between the 
frequency response of the system and the frequency response of 
its estimate. 
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Since we try to fit a finite impulse response model to a 
recursive system (with impulse response of infinite length) , 
we cannot pretend to estimate the frequency response at all 
frequencies in an effective way. Therefore, we input test 
signals at finite number of frequencies and we look at the 
estimated spectrum at these frequencies only. 

As it will be seen in the sequel, the results are 
satisfactory even in the presence of added disturbances in the 
measurements. The recursive use of blocks of data has the 
effect of smoothing the effect of disturbances in the 
estimates. 

D. SIMULATION RESULTS 

In this thesis, we try inputs with different frequencies 
to test the identification algorithm. Also, we test the 
behavior when the input is the combination of several 
different frequencies to examine the effectiveness of the 
technique. In the next paragraph, we discuss some results 
obtained by running the program. 

1. Original System 

Figure 3.1 shows the magnitude and phase of frequency 
response of the original system. It shows the low pass nature 
of the system. 
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2. Single Frequency Input at Arbitrary Frequency without 
Noise 

Figure 3.2 and 3.3 are in the same group in which we 
tried a constant input with magnitude 10. As predicted, 
figure 3.2 portrays the error EPS between the output and 
predicted output which decreases and approaches zero. 
Furthermore, in figure 3.3, the frequency response of the 
estimated model at the input frequency is shown to be 
identical to the one of the original system. 

In a second set of experiments, we tried several 
sinusoidal inputs with different frequencies. First of all, 
we used an input signal with a digital frequency tt/2 and 
obtain the results shown in figures 3.4 and 3.5. And then, we 
tried another input with frequency tt/8 and get the results 
shown in figures 3.6 and 3.7. In like manner, the error 
approaches to zero quickly and the frequency response of the 
estimated model at the input frequency is the same as the one 
of the orignal system. 

3. Multiple Frequency Inputs without Noise 

In a third set of experiments, we tried the 
combination of four frequency inputs at arbitrary frequencies 
to obtain the results shown in figures 3.8, 3.9, and 3.10. 

The output in the time domain is shown in figure 3.8. In this 
case, the error between true and predicted outputs approaches 
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zero quickly and the frequency responses at the input 
frequencies are still identical. 

4. Single Frequency Input with Noise 

In this section we show results of similar experiments 
with random measurement noise. The noise considered is 
Gaussian white zero mean, i.e., independently and indentically 
distributed. Its standard deviation is set at about one third 
of the output magnitude . 

As before we first excite the system with a constant 
input, and we add measurement noise to the output of the 
system. The results are shown in figures 3.11, 3.12, and 
3.13. From figure 3.12, even though we find that the error 
does not decay to zero due to the presence of noise, the 
frequency response of the estimated model is sufficiently 
close to the one of the original system. 

In a second set of experiments, we tried sinusoidal 
inputs again with noisy measurements. There are four inputs 
with different digital frequencies (llir/16 ,ir/2 , 7r/4, and n/16) 
used and the results are shown in figures (3.14-3.22) 
respectively. Then, by inspecting the figures, the error 
still does not go to zero due to the presence of noise, but 
the estimated frequency response is sufficiently close to the 
one of the original system. 
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5. Multiple Frequency Inputs with Noise 

In like manner, we input signals with multiple 
frequencies to test the alogrithm. In this case, we used four 
different groups of inputs. The first one is constant and 
sinusoidal (frequency n/2) input. The 2nd one is a constant 
and sinusoidal (frequency n/4) input. The 3rd and 4th are 
sinusoidal with two frequencies (7r/2,7r/4) and (7r/4 ,7t/8) 
respectively. The results are shown in figures 3.23 to 3.31. 
Then, inspection of these figures reveals that the error in 
each case is bounded and the frequency responses of the 
estimated model are sufficiently close to the one of the 
original system. 
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FREQUENCY RESPONSE OF TUB ORIGINAL SYSTEM H(Z) 
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THETA 

Figure 3.1 Frequency Response of the Original System 



x 10 3 EPS (when u(n)= 10 y(n)=yp(n) 1-0* rand ) 
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H(Z) & ESTIMATED H(Z) AT THE INPUT FREQUENCY 
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EPS (when u(n)=20cos(pi*n/2) y(n)=yp(n) l O*rand ) 
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Figure 3.4 Error Between Original and Predicted Output 



H(Z) & ESTIMATED I I(Z) AT THE INPUT FREQUENCY 
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EPS (when u(n)=40cos(pi*n/8) y(n)=yp(n) HO’rand ) 
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Figure 3.6 Error Between Original and Predicted Output 
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YP(n) ORIGINAL OUTPUT WITHOUT NOISE 
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Figure 3.8 The Original Output without Noise 



EPS (when u(n)= 10+20cos(pi*n/2)+30cos(piW4)+40cos(pi*n/8) y(n)=yp(n) ) 
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Figure 3.9 Error Between Original and Predicted Output 



H(Z) & ESTIMATED H(Z) AT THE INPUT FREQUENCY 
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Figure 3.10 The Original and Estimated Frequency Response 



YP(n) ORIGINAL OUTPUT WITHOUT NOISE 
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Figure 3.11 The Original Output with Noise and without Noise 



EPS (when u(n)=10 y(n)= yp(n)+ 25*rand ) 
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Figure 3.12 Error Between Original and Predicted Output 



H(Z) & ESTIMATED II(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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YP(n) ORIGINAL OUTPUT WITHOUT NOISE 
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Figure 3.14 The Original output with Noise and without Noise 



EPS (when u(n)=20cos(lPpi*n/16) y(n)=yp(n)t-35*rand ) 
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Figure 3.15 Error Between Original and Predicted Output 



1I(Z) & ESTIMATED 1I(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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THETA ( PI*K/(N/2) ) 

Figure 3.16 The Original and Estimated Frequency Respons 



EPS (when u(n)=20cos(pi*n/2) y(n)= yp(n)+35*rand ) 
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Figure 3.17 Error Between Original and Predicted Output 



H(Z) & ESTIMATED II(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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EPS (when u(n)=30cos(pi*n/4) y(n)=yp(n) P60*rand ) 
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Figure 3.19 Error Between Original and Predicted Output 



H(Z) & ESTIMATED H(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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EPS (when u(n)= 50cos(pi*n/16) y(n)=yp(n)+ 120* rand ) 
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Figure 3.21 Error Between Original and Predicted Output 



I I(Z) & ESTIMATED 11(Z) AT TI IE INPUT FREQUENCY (WITH NOISE) 
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EPS (when u(n)=10+20cos(pi*n/2) y(n)=yp(n)+50*rand ) 
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Figure 3.23 Error Between Original and Predicted Output 



I'I(Z) & ESTIMATED 1I(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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TI IETA ( PI*K/(N/2) ) 



EPS (when u(n)=10l-30cos(pi*n/4) y(n)= yp(n) l-90'rand ) 
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Figure 3.25 Error Between Original and Predicted Output 



1I(Z) & ESTIMATED II(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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THETA ( PI*K/(N/2) ) 

Figure 3.26 The Original and Estimated Frequency Response 



EPS (when u(n)=20cos(pi t n/2)+30cos(pi*n/4) y(n)=yp(n)H-90*rand ) 
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Figure 3.27 Error Between Original and Predicted Output 



H(Z) & ESTIMATED I1(Z) AT THE INPUT FREQUENCY (WITH NOISE) 
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YP(n) ORIGINAL OUTPUT WITHOUT NOISE 
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Figure 3.29 The Original Output with Noise and without Noise 



EPS (when n(n)=30cos(pi t n/4)+40cos(pi*n/8) y(n)=yp(n)t- 140* rand ) 
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Figure 3.30 Error Between Original and Predicted Output 



I I(Z) &. ESTIMATED II(Z) AT TI IE INPUT FREQUENCY (WITH NOISE) 
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IV. CONCLUSION 



In this thesis we introduced an approach to the 
identification of linear models based on the frequency domain 
formulation. The recursive nature and the use of FFT 
techniques make the approach attractive for on-line 
implementation . 

A particular aspect we addressed is the identification of 
the frequency response of the system by approximating the 
impulse response with a finite sequence. By use of a 
simulated example we have explored the convergence properties 
of the algorithm, and its robustness in the presence of 
measurement noise. 

Several issues remain to be investigated, mainly the 
advantages (if any) of this approach in comparison with more 
conventional time domain estimation techniques. Although this 
is subject of future research, as a conjecture we can say that 
the recursive frequency domain approach might be more robust 
than the conventional time domain in the cases when the input 
signal has energy concentrated around a few frequencies, such 
as the case of periodic excitation. In this case the 
processing gain of the DFT should give a better performance in 
the presence of measurement noise. 
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APPENDIX A : MAINPROGRAM 



. 9. q. O. O. O. 



.9.9.9.9.9.9.5.S.9.9-9-9.9-9-S-5-9-3.9-9.9.9-9.9.9-9-9-9-5-9-9-5-9.9.9.9.9.9.9.9.9.9.9.9.9.4.9.9.9.0. 



-o 

% THESIS MAIN PROGRAM File name : THESIS4 . M 

% 



% Advisor : Roberto Cristi 
% Student : Chao, Chi-Shun 
% Date : Jan, 1, 1991 



'd'b'b'o'o'o'o'o'o'o'o'o'o'o'o'o'o^'o'o'o'o'o'o'o'o^'^ 



- 2* 2- S' 2- 5- 2- 5- 5- 

'b'b'b'b'b'o'b'6 



j'o'b'O'&'o'b'b'o'b'b'o 



% Set some constant number 



clg 

clear 

tfinal = 80; 
dt = 0.1; 
kmax = tfinal/dt; 
rand ( 1 normal 1 ) ; 
rand_mag = 35; 

M = 20; 

N = 32; 

N2= N/2 ; 

% Set the initial condition 
yp(l)=0; 

YP ( 2 ) =0 ; 

y ( i) =° ; 

y ( 2 ) =o ,* 
u ( 1) =0 ; 
u ( 2 ) =0 ; 



% The max. time to run 
% Sampling time 
% The max. number of signals 
% Set the type of noise 
% Set the gain of noise 
% The # to calculate the error 
% The # of signals in a block 



% Produce the input and coresponding output 
% at different frequencies 
for n=3 : kmax 

%u(n)= 10; 

u(n)= 20*cos (pi*n*5/16) ; 
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o\° o\° o\° o\o o\o o\° o\o o\° 



%u(n)= 30*cos(pi*n/4) ; 

%u(n)= 40*cos (pi*n/8 ) ; 

%u(n)= 50*cos (pi*n/16) ; 

%u(n)= 60*cos (pi*n/32) ; 

%u(n)= 10+30*cos(pi*n/4) ; 

%u(n)= 40*cos (pi*n/8 ) +30*cos (pi*n/4 ) ; 

%u(n)= 10+20*cos(pi*n/2)+30*cos(pi*n/4)+40*cos(pi*n/8) ; 
%u(n)= 10+20*cos (pi*n/2 ) +30*cos (pi*n/8 ) ... 

+40*cos (pi*n/8+50*cos (pi*n/16) ; 

yp (n) =0 . 3*yp (n-l)+5.0*u(n-l) ; 
y(n) =yp ( n ) + r and_mag * rand ; 



end 



% Produce the figures of input and output 
plot (u) ; 

title ( 'INPUT U(n) ' ) ; 

xlabel ( 'n' ) ; ylabel( 'MAGNITUDE') ; 

grid; 

pause ; 

clg 

subplot(211) ,plot(yp) ; 

title ( ' YP (n) ORIGINAL OUTPUT WITHOUT NOISE '); 

xlabel ( ' n ' ) ; ylabel ( ' MAGNITUDE ' ) ; 

grid; 

subplot(212) ,plot(y) ; 

title ( ' Y (n) OUTPUT WITH NOISE U (n) =20cos ( 5*pi*n/ 16 ) ... 

y(n)=yp(n) ' ) ; 

xlabel ( ' n ' ) ; ylabel ( ' MAGNITUDE ' ) ; 

grid; 

%meta thesis44 
pause ; 
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% Produce the frequency response of the original system 

% by using the "dbode" function 

num=[5] ; 

den=[l,-0.3] ; 

w=linspace(0,pi, 17) ; 

[mag, phase ] =dbode (num, den, w) ; 

clg 

subplot(211) ,plot(w,mag) ; 

title ( 'FREQUENCY RESPONSE OF THE ORIGINAL SYSTEM H(Z)') 
xlabel ( 'THETA' ) ; 
ylabel ( 'MAGNITUDE' ) ; 
grid; 

subplot(212) , plot (w, phase) ; 

title (' FREQUENCY RESPONSE OF THE ORIGINAL SYSTEM H(Z)') 

xlabel ( 'THETA' ) ; 

ylabel ( 'PHASE (DEGREE) ') ; 

grid; 

pause; 



% Create the matrix H 
v= [ zeros (1,N2) ,ones(l,N2) ] ; 
h=diag(v) ; 

F=f ft (eye (N) ) ; 

H=F*h*inv (F) ; 



% Set the zero matrix 
THETAK = zeros (N , 1) ; 
P0=10*eye (N) ; 

P=PO ; 
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% Use the FFT and RLS (Recursive Least Square) techniques 
% to estimate frequency response of the model. It is an 
% interactive and iterative loop 
for k = 2 : M 

ykO = y ( (k-1) *N2+1 : k*N2 )'; 

uk = u ( ( k-2 ) *N2+1 : k*N2 )'; 

yk=[zeros (N2 , 1) ;ykO] ; 

YK = fft( yk ) ; 

UK = fft ( uk ) ; 

EPS (k) =0.0; 

for i = 1 ; N 

PHI=diag (H ( i , : ) ) *UK ; 

[ THETAK , P ] =r 1 s ( THETAK , P , PHI , YK ( i ) ) ; 

EPS (k) =EPS (k) + ( YK (i) -conj (PHI') *THETAK) '* 
(YK(i) -conj (PHI ' ) *THETAK) ; 

end 



THETA ( : , k) = THETAK ; 



end 

EPS = sqrt (EPS) /N; 

% Calculate the corresponding horizontal axis 
for k=0 ; N/2 

wl (k+1) =2*pi*k/N ; 

end 
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% Produce the frequency response figures of estimated model 
clg 

subplot (2 11) ,plot (wl,abs(THETA(l: 17, M) ) ) ; 

title ( ' EST . H (Z ) ' ) ; 

xlabel( 'THETA ( PI*K/(N/2) )'); 

ylabel ( 'MAGNITUDE ' ) ; 

grid; 

subplot (2 12) , plot (wl, angle (THETA(1: 17 ,M) ) ) ; 

title ( 'EST. H ( Z ) ' ) ; 

xlabel ( 'THETA ( PI*K/(N/2) )'); 

ylabel ( 'PHASE (RADIAN)'); 

grid; 

pause ; 



% Produce the error figures 
clg 

displot (EPS ) ; 

title('EPS (when u(n) =20cos (5*pi*n/16) 
y (n) =yp (n) +35*rand )'); 
xlabel ( 'K') ; 
ylabel ( 'MAGNITUDE ' ) ; 
grid; 

%meta thesis44 
pause ; 



% Produce the frequency response of the estimated model 

% in a continuous pattern. 

clg 

plot (w, mag , wl , abs (THETA (1:17, M))); 
title ( ' H (Z ) & EST. H(Z) '); 

xlabel ( 'THETA ( PI*K/(N/2) )'); 

ylabel ( ' MAGNITUDE ' ) ; 
grid; 
pause ; 
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% Store the value of the response of the estimated model 
% at the input frequency. 

THETA_F = zeros (17,1) 



%THETA_F (1,1) 
%THETA_F (2,1) 
%THETA_F ( 3 , 1) 
%THETA_F (5,1) 
%THETA_F (9,1) 
THETA_F (6,1) 



THETA (1,M) 
= THETA ( 2, M) 
= THETA ( 3, M) 
= THETA ( 5, M) 
= THETA ( 9, M) 
THETA ( 6, M) ; 



% Produce the frequency response of the original system and 
% the response of estimated model at the input frequency 
% in the same figure, 
clg 

plot (w, mag) ; 

title ( ' H ( Z ) & ESTIMATED H(Z) AT THE INPUT FREQUENCY ... 

(WITH NOISE) ' ) ; 

xlabel ( 'THETA ( PI*K/(N/2) )'); 
y label ( ' MAGNITUDE ' ) ; 
grid; 
hold on 

displot (wl , abs (THETA_F) ) ; 
hold off 
%meta thesis44 
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o\o o\o o\o o\o o\o o\o o\o o\o <U o\o o\° o\° o\o o\° o\o o\o o\o o\o o\° o\° o\o o\° o\o o\o 



APPENDIX B 



FUNCTION TYPE SUBROUTINE RLS 

I 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 



SUBROUTINE TYPE 



FUNCTION 



File name : RLS.M 



Advisor 

Student 

Date 



Roberto Cristi 
Chao, Chi-Shun 
Jan, 1, 1991 



'o 'o '6'b 'o'o'o'o'o'o^'b'o'b'o'b^'o^^'o'o'o'b^'o'o'^'o'b'^'o'o'6'^'b'o'b'o^'S'b'o'o'b'o'o'b'b'b'o'b'o'o'o'o'b'o 



unction [x, p] =rls (x, p, phi 



Y) 



'o'S'o'o'b'o'S'o'b'b^'o^^'o^'o'o'^^'o'o'o'b'b^^'o^'o'o'o'o'b^'o'o'o'S'o'o'o^^'o'o'o'o^'o^'o^'b'o'o'S'b 



It updates the estimate x and its covariance p 
by standard recursive least squares. 

data: x = column vector (input, output) 

p = square matrix, positive definite ( input , output) 

y = scalar (input) 

phi = column vector (input) 

NOTE: Due to numerical approximations it might be 

a good idea to check the matrix p for positive 
definiteness . 



%%%%%%%%%%%%%%%%%%: 

oooooooooooooooooo 



oooooooo A 






'b'o'o 'b ' o 



fact = 1.0 + phi'*p*phi; 

x = x + p*conj (phi) * (y-conj (phi ') *x) /fact ; 
p = p - p*conj (phi) *conj (phi' ) *p/fact ; 

% end rls 
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